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The nuclear velocity perturbation theory (NVPT) for vibrational circular dichroism (VCD) is de¬ 
rived from the exact factorization of the electron-nuclear wave function. This new formalism offers 
an exact starting point to include correction terms to the Born-Oppenheimer (BO) form of the molec¬ 
ular wave function, similarly to the complete-adiabatic approximation. The corrections depend on a 
small parameter that, in a classical treatment of the nuclei, is identified as the nuclear velocity. Apart 
from proposing a rigorous basis for the NVPT, we show that the rotational strengths, related to the 
intensity of the VCD signal, contain a new contribution beyond-BO that can be evaluated with the 
NVPT and that only arises when the exact factorization approach is employed. Numerical results are 
presented for chiral and non-chiral systems to test the validity of the approach. 


I. INTRODUCTION 

Vibrational circular dichroism (VCD) fll-lTl in 
molecules refers to the difference in absorption of 
left and right circularly polarized light in the infrared 
region of the electromagnetic spectrum. In contrast to 
circular dichroism, that originates in electronic transi¬ 
tions, VCD is the difference in interaction of a molecule 
with radiation of opposite circular polarizations when 
it undergoes vibrational transitions. Experimentally, 
VCD is employed to probe the absolute configuration 
of chiral molecules in solution and provides detailed 
structural information, thus being a very sensitive form 
of vibrational spectroscopy. 

From the theoretical point of view |43-j2jj, the Born- 
Oppenheimer (BO) l22ll . or adiabatic, treatment of the 
coupled motion of electrons and nuclei in molecular sys¬ 
tems is inadequate for predicting VCD. Since the inten¬ 
sity of the VCD signal is proportional to the rotational 
strength for a transition between two vibrational states, 
the calculation of the electric current and of the mag¬ 
netic dipole moment (and of their scalar product) is re¬ 
quired. The electric current and the magnetic dipole mo¬ 
ment contain both electronic and nuclear contributions, 
but when the BO approximation is employed, the elec¬ 
tronic contributions identically vanish. This is due to 
the fact that the ground-state electronic wave function 
is real for a non-degenerate adiabatic state and therefore 
the expectation values of the purely imaginary (Hermi- 
tian) electric current |23ll28l and magnetic dipole mo¬ 
ment operators vanish ||l8j. Therefore, VCD appears a 
fundamentally non-adiabatic (beyond-BO) process, thus 
requiring a theoretical approach able to explicitly treat 
the dynamical coupling between electronic and nuclear 


degrees freedom in molecules. 

A practical question l29i arises at this point, as to 
whether such coupling can be accounted for within 
a standard ab-initio molecular dynamics formulation. 
Among the most successful ideas are in fact those re¬ 
sorting to the treatment of beyond-BO effects as a per¬ 
turbation to the BO problem, numerically less expen¬ 
sive than a full non-adiabatic calculation but indeed not 
consistent if strong non-adiabatic effects are expected, 
e.g. in the presence of conical intersections. Exam¬ 
ples are the approaches proposed by Nafie lfT9l , employ¬ 
ing the complete-adiabatic expression of the electron- 
nuclear wave function, and by Stephens If20ll . introduc¬ 
ing the magnetic field perturbation theory. These meth¬ 
ods allow to overcome the problems encountered in 
the BO calculation of VCD, while exploiting the ad¬ 
vantages of the BO formalism like the product form of 
the electron-nuclear wave function. Recently, VCD has 
been calculated by developing and implementing a nu¬ 
clear velocity perturbation theory (NVPT) f30l based on 
the complete-adiabatic approach of Nafie lfl9ll . In this 
formulation, non-adiabatic corrections to the electronic 
adiabatic ground-state are perturbatively taken into ac¬ 
count and are induced by a "small" nuclear velocity. 

In this paper we propose a novel approach to NVPT, 
based on the exact factorization of the electron-nuclear 
wave function Enna. The advantage of this formula¬ 
tion comes from using a product form, like in the BO 
approximation, of the wave function, that is not the 
result of an approximation but an exact starting point. 
The electron-nuclear wave function is written as a sin¬ 
gle product of an electronic many-body factor, para¬ 
metrically depending on the nuclear positions, and a 
nuclear wave function. The latter can be interpreted 


as a proper nuclear wave function since it leads to 
the exact nuclear density and current-density. More¬ 
over, when the product form is inserted into the time- 
dependent Schrodinger equation (TDSE), two coupled 
equations for the components of the full wave function 
are derived, with the nuclear equation being a TDSE 
where a time-dependent vector potential and a time- 
dependent scalar potential (or time-dependent poten¬ 
tial energy surface, TDPES) 133H371 represent the effect 
of the electrons on the nuclei beyond-BO. Therefore, in 
this context, the electronic equation generates the proper 
evolution expected when the coupling between elec¬ 
trons and nuclei is fully accounted for and it allows to 
recover the BO equation in a certain limit, as will be dis¬ 
cussed in the paper. 


Two major results will be reported: (i) NVPT l30ll will 
be rigorously derived, using as starting point the exact 
electronic equation from the factorization rather than 
the complete-adiabatic approach 119), and (ii) correc¬ 
tion terms to the "standard" expression of the rotational 
strength will naturally appear in the new formulation, 
due to the presence of the time-dependent vector poten¬ 
tial of the theory. Throughout the paper, we will adopt 
a time-dependent picture, as this is crucial to introduce 
the concept of nuclear velocity and, thus, to make the 
connection with NVPT. In such a time-dependent pic¬ 
ture we will have access to the instantaneous expecta¬ 
tion values of the electric current and of the magnetic 
dipole moment. The corrections to those expectation 
values, and therefore to the rotational strength, can be 
derived also in a static picture referring to the time- 
independent formulation lf38l of the factorization, but 
the direct link to NVPT would then be missing. 


IIA 


The paper is organized as follows. In Section 
we review the linear response theory approach to VCD, 
showing the connection between rotational strength and 
intensity of the spectrum. In Section [II B| we recall the 
exact factorization formalism. In Section [III] we focus 
on the electronic equation from the exact factorization, 
showing how to recover the BO limit and introducing 
non-adiabatic effects as a perturbation to the adiabatic 
framework. The perturbation parameter is identified as 
the nuclear velocity, exactly as in NVPT, if the classical 
limit is considered. However, here we have access to 
the quantum electronic evolution equation, thus the per¬ 
turbation parameter has a more general meaning since 
we are not restricted to a classical treatment of the nu- 



d 2 -oxirane, a chiral system, in comparison to oxirane, a 
non-chiral molecule. We also report the comparison be¬ 
tween the NVPT approach and the more standard mag¬ 
netic field perturbation theory l20l in Section|V B[ for (S)- 


d 2 -oxirane, (R)-propylene-oxide and (R)-fluoro-oxirane. 
Our conclusions are stated in Section IVTI 


II. THEORETICAL BACKGROUND 
A. Vibrational circular dichroism 

Vibrational spectroscopy probes the coupling of the 
nuclear degrees of freedom to applied electro-magnetic 
fields. On the macroscopic level, the absorption process 
is described phenomenologically by the Beer-Lambert 
law [|39l , where the material specific attenuation of the 
radiation per unit length is accounted for by the molar 
absorption coefficient e. Microscopically, and within the 
linear response regime, the energy dissipated in the in¬ 
teraction between the medium and the radiation is ex¬ 
pressed in terms of the observable that couples to the 
external field. In case of radiation in the infrared spec¬ 
tral range, the multipole approximation and the long 
wavelength limit can be applied l39l [40l to determine 
such coupling. The microscopic and the macroscopic 
perspectives can be connected in the framework of lin¬ 
ear response theory ETI . In the Heisenberg formu¬ 
lation, the frequency dependent absorption coefficient 
takes the form of the power spectrum of the dipole auto¬ 
correlation 14211431. 

The specific feature of VCD is the different interaction 
of chiral systems with polarized light. Linearly polar¬ 
ized light encounters optical rotatory dispersion while 
circularly polarized light encounters different absorp¬ 
tions for the different handednesses of the radiation, the 
vibrational circular dichroism (VCD). Formally, this is 
accounted for by the dependence of the refractive index 
of a chiral system on the handedness of the radiation. 
While this effect is not relevant for the mean infrared ab¬ 
sorption, the difference absorption gives rise to the VCD 
signal. 

For the calculation of the absorption coefficient e(u>) 
and of the difference absorption A e(ui) 0, a common 
approach in the literature is to invoke the double har¬ 
monic approximation for nuclear motion and dipole 
moment. This leads to the expressions 

and 

o_ 3 

Ae( “ ) = 4 3^o,( M ) N fl ^-^). (2) 

The dipole strength D *. and rotational strength Iij ; of the 
vibrational mode k, with frequency u>k, are evaluated as 


D„ = »<?> 
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respectively, with the time derivative of the dipole mo¬ 
ment fi, namely the current, and the magnetic dipole 
moment rh. In Eqs. (T) and (5), V indicates the volume 
occupied by the system, h = 2nh is the Planck constant, 
c is the speed of light, and n(iS) is the refractive index of 
the medium. Normal modes will be indicate throughout 
the paper as q, with velocities q. 

The linear variations of the expectation values (over 
the instantaneous state of the system) of the current and 
of the magnetic dipole moment with respect to (w.r.t.) 
the mode q k around their equilibrium values are calcu¬ 
lated from the total (electronic and nuclear) atomic po¬ 
lar tensor V v (APT) and atomic axial tensor M. 1 ' (AAT). 
The APT and AAT have electronic and nuclear contribu¬ 
tions 1311301, namely 


djdfi) _ rrjl! _ CV 

dRz 0/3 “ “ 0 


■Ml 


Gift 


djrhp) 

dK 


= M 


al3 


%a0 + Kpi 
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( 6 ) 


with electronic parts £ and X and nuclear parts M and 
J. Here, the indices a and /3 are used for the Cartesian 
coordinates, while v labels the nuclei. The dipole and 
rotational strengths are related via the chain rule to the 
vibrational nuclear displacement vector S v ak which de¬ 
scribes the displacement of nucleus v in direction a due 
to the fc-th normal mode q k , 


™ = dR £ 
“ fe dq k 


q=0 


dK 

dqk 


q—0 


(7) 


B. Exact factorization of the electron-nuclear wave 
function 


can be exactly factorized to the product 

T(r,R,f) = <t>R,(r, f)x(R, t) 


( 11 ) 


where 


/ 


dr 


l$R(r,t)| 2 


1 V R, t. 


( 12 ) 


Here, y(R. t) is the nuclear wave function and $ R (r, t) 
is the electronic wave function which parametrically de¬ 
pends on the nuclear positions and satisfies the partial 
normalization condition (PNC) expressed in Eq. (12) . 
The PNC guarantees the interpretation of |\(R, f)| cis 
the probability of finding the nuclear configuration R at 
time t, and of |$ R (r,f)| 2 itself as the conditional prob¬ 
ability of finding the electronic configuration r at time 
t, given the nuclear configuration R. Further, the PNC 
makes the factorization 111' unique up to within a (R, in¬ 
dependent gauge transformation. 


X(R, t) x(R, t) = e fi^ R ’*)x(R,t) 
$R(rd)-t ( i , R(r,t) = eR^ R ^)$ R (r,i). 


(13) 


where 0(R, t) is some real function of the nuclear coor¬ 
dinates and time. 

The stationary variations j44l of the quantum me¬ 
chanical action w.r.t. c l > R (r,i) and ;yfIT., t) lead to the 
equations of motion 

[Helix, R) - e(R, t )) $ R (r, t) = ihd t $ R (r, t) (14) 
H„ (R, t)x(R, t) = ihdtxi R, t), (15) 

where the PNC is inserted by means of Lagrange multi¬ 
pliers MM- Here, the electronic and nuclear Hamilto¬ 
nians are defined as 


Kii r, R) = Hbo( r, R) + U™ up [$r, X ] (16) 


The non-relativistic Hamiltonian describing a system 
of interacting electrons and nuclei, in the absence of a 
time-dependent external field, is 


H = T n + H B o, 


( 8 ) 


where T n is the nuclear kinetic energy operator and 


and 


H n (R,t) 


[-iftMv + A„(R, t )] 2 

^-2 W v - +e(R ’ 0 

v—1 


(17) 


respectively, with the "electron-nuclear coupling opera¬ 
tor" 


H BO {r , R) = T e (r) + W ee ( r) + V en (r, R) + W nn ( R) (9) 

is the standard BO electronic Hamiltonian, with elec¬ 
tronic kinetic energy T e (r), and with potentials W ee ( r) 
for electron-electron, W nn ( R) for nucleus-nucleus, and 
14 n (r,R) for electron-nucleus interaction. The symbols 
r and R are used to collectively indicate the coordinates 
of N e electrons and N n nuclei, respectively. 

It has been proved |3ll [32 1 that the full time- 
dependent electron-nuclear wave function 'T(r,R, t) 
that is the solution of the TDSE, 


CT*r,x] = £m; l 


Nn 1 r[-av„-A„(R,o ] 2 


(18) 


X 


+ A„(R, t) ) (-ifiV v - A„(R,i)) 


The time-dependent potentials are the TDPES, e(R, £), 
implicitly defined by Eq. (14) as 

e(R,t) = ($ R (f)| H bo + Ue° up - ihd t |$ R (f)) r , (19) 
and the vector potential. A,, (R. /:), defined as 


R, t) = ihd t ^(r, R, t), (10) 


A, y (R, t) = ($ R (t)| -ifiV, $ R (i)) r • (20) 
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The symbol ( • ) r indicates an integration over electronic 
coordinates only. Under the gauge transformation (13) , 
the scalar potential and the vector potential transform as 


e(R,i) =e(R,t)+&0(R,i) (21) 

A„(R, t) = A„(R, f) + V v 0(R, t). (22) 


In Eqs. 114 1 and 115 , Ug°“ p [$ R , \], e(R, t), and A„ (R, t) 
are responsible for the coupling between electrons and 
nuclei in a formally exact way It is worth noting that 
the electron-nuclear coupling operator, U™“ P [$ R ,\] in 
the electronic equation (14} depends on the nuclear 
wave function and acts on the parametric dependence 
of $ R (r,t) as a differential operator. This "pseudo¬ 
operator" includes the coupling to the nuclear sub¬ 
system beyond the parametric dependence in the BO 
Hamiltonian 7T R o(r, R). 

The nuclear equation (15} has the particularly appeal¬ 
ing form of a Schrodinger equation that contains the 
TDPES (19} and the vector potential (20} governing nu¬ 
clear dynamics and yielding the nuclear wave func¬ 
tion. The scalar and vector potentials are uniquely de¬ 
termined up to within a gauge transformation, given by 
Eqs. (21} an d (22} . As expected, the nuclear Hamilto- 
is form-invariant under such transfor¬ 


man in 


mations. x(R, t) is interpreted as the nuclear wave func¬ 
tion since it leads to an A'-bod y nuclear density and an 
A-body current-density which reproduce the true nu¬ 
clear A-body density and current-density l32l obtained 
from the full wave function ^(r, R, t). The uniqueness 
of e(R, t) and A„(R, t) can be straightforwardly proved 
by following the steps of the current-density version [47J 
of the Runge-Gross theorem [48], or by referring to the 
theorems proved in Ref. [31]. 


III. NUCLEAR VELOCITY PERTURBATION THEORY 

Before showing the derivation of the velocity- 
dependent corrections to the BO wave function within 
the exact factorization approach, let us present a pro¬ 
cedure that allows us to recover the BO limit of the 
electronic equation (14} . Suppose first that the electron- 
nuclear wave function is given as a BO product 


The first term on the right hand side (RHS), containing 
the Laplacian 15114531 w.r.t. nuclear coordinates, will be 
neglected from now on. It can be shown, as reported 
in Refs. 15444561 . that this term contributes with second- 
order non-adiabatic couplings to the electronic equation, 
but being explicitly 0(M~ 1 ) its effect can be neglected if 
compared to the remaining (and leading) term. Follow¬ 
ing again Refs. 15444561 1. the term that depends on x can 
be approximated to zero-th order in an //-expansion l57l 
of the nuclear wave function as the classical nuclear ve¬ 
locity, namely 

M v x(R, t) M v 


We have invoked here the classical limit in order to di¬ 
rectly relate our results to the NVPT Il30l and to justify 
the condition of "small nuclear velocity" that allows a 
treatment of effects beyond-BO within perturbation the¬ 
ory. The procedure, however, does not rely on the clas¬ 
sical limit and the "small" perturbation parameter will 
be denoted as 


Ai,(R, t) 


1 -fhV„x(R. t) 
M v x(R ,t) 


(26) 


Eq. (26} contains the variations in space of the phase and 
of the modulus of the nuclear wave function 15B1 , and 
when both variations are "small" then the approach con¬ 
sidered here can be applied. We have justified the for¬ 
mer hypothesis (small variations of the phase) by em¬ 
ploying the classical approximation and we are now as¬ 
suming valid also the latter (small variations of the mod¬ 
ulus). 

The electronic Hamiltonian from Eq. (16} becomes 


JV„ 


H el ( r, R) = H bo + ^ A„(R, t) • (-%KV„) (27) 


and the TDPES reads 


N n 


e(R, t) = H B o + ^ A„(R, t) • (-ihX7 u ) 

V=l 

= 4°o(R)> 


( 0 ) \ 
^ ; r 


(28) 


4 (r, R, t) = $ R (r, t)\ (R, t) = (r) x (R, t). (23) 


Here, <p R ! (r) indicates the real and not degenerate BO 
ground-state. Using Eq. (20} , it follows that the vector 
potential vanishes identically, A„(R, t) = 0 |[49| . This 
can be interpreted as a choice of gauge [501. With this as¬ 
sumption, the electron-nuclear coupling operator from 
Eq. (18} becomes 


N n 

U C e°n UP \^X\=Y, 

v =1 


-h 2 vl 

2 M v 


-ffiv^x( R u) 

M v \{ R, t) 


(—ihV„) ■ 

(24) 


i.e. only the Hbo term survives, since the second term 
does not contribute to the TDPES. Notice that here the 
term (4 > R (t)|— ihdt\^n(t)) r identically vanishes, because 
the electronic wave function is the time-independent BO 
wave function. In order to recover from Eq. (27} the 
electronic equation within the BO approximation, one 
should impose A„(R, t) = 0, or similarly R v (t) = 0 V v 
as the electronic equation in BO is solved for fixed nuclei 
(meaning that their velocity is zero). 

To summarize, in order to construct the Hamiltonian 
in Eq. (27} (i) we treat the nuclei classically, thus we 
considerthe nuclear wave function up to within O(h 0 ) 
terms, (ii) we derive corrections to the BO Hamiltonian 


4 











that are proportional to the nuclear velocity, thus recov¬ 
ering the BO electronic equation if the nuclear velocity 
is zero (condition of fixed nuclei), (iii) we "relax" the 
hypothesis of classical nuclei by introducing A„(R, t) as 
the perturbation parameter. 

Combining Eqs. 1 27] and [28] will provide the elec¬ 
tronic equation within the new formulation of NVPT. 
In contrast to the formulation based on the complete- 
adiabatic approach |[19l , the perturbative scheme pre¬ 
sented here directly applies to the electronic equation 
rather than to the full TDSE. Using perturbation the¬ 
ory [30], where Hbo is the unperturbed Hamiltonian 
and the second term on the RHS of Eq. [27] is the per¬ 
turbation, we find the solutions of the equation 


dependence on Aj,(R, t), and <PRra( r ) is orthogonal to 

yp 1 ' (r). This last property can be interpreted as a choice 
of gauge. For instance, by imposing the condition that 
(<p^|<PR,(t)) r is real VR,t, which is allowed as gauge 
condition, we imply the orthogonality of yp| :, J/(t (r) and 

<Pp^( r), namely (<PR^|<PR^„ a ) r = 0. It easy to prove that 
the PNC remains valid up to within 0( A"), using the 
orthogonality of ^'(r) and 

The first-order approximation to the TDPES is 

e (1) (R, t) = (R) - * Y, °(M R , *)) (33) 


N n 

v=\ 


p R (r,t) = 0, 


(29) 


as 


<PR.(r, t) = <p R ’ (r) 


E 

e^O 


( e ) 




¥>R ( r ) 

(30) 


e BO (R) — e BO (R) 


up to within linear-order in the perturbation, with the 
index v running over the N n nuclei and with a running 
over the three Cartesian coordinates. The symbol 9" is 
used to indicate a spatial derivative along the a direc¬ 
tion of the position of the v-\h nucleus and e labels the 
(unperturbed) adiabatic excited states. The TDPES, up 
to within first-order terms, is labeled b'k It is worth not¬ 
ing that in writing Eq. [29] , we are discarding the varia¬ 
tions in time of the first-order correction to the BO wave 
function, adopting a previously assumed l30l hypoth¬ 
esis that these variations are smaller than the perturba¬ 
tion itself, thus negligible at the given order. We re-write 
Eq. (|30j as 

<m(r,t) = p R (r) +E* A “( R ’^^Ea( r )> (3!) 

V,OL 

introducing the definition of the first-order perturbation 
to the BO ground-state 

(1) / \ V - '' d e o va (R-) (e) / 

< 32 > 

Here d,, 0 „ o (R) is the a-th Cartesian component of the 
non-adiabatic coupling vector, corresponding to the v- 
th nucleus, between the unperturbed ground-state and 
the excited state e, whereas the frequency w e0 (R) is the 
the energy difference (divided by h) between the excited 
(e) and the ground (0) states. When the adiabatic states 
are real, Eq. [32] is real as well and the second term in 
Eq. j3T) is purely imaginary. Moreover, the correction 
term in Eq. [3l] depends on time only implicitly, via its 


but the second term on the RHS is identically zero, as 
can be proved by either inserting Eq. ( |31~] in the defini¬ 
tion of the TDPES given in Eq. ( [T9] ) orby considering 
the fact that e^ 1 ) (R, t) must be real while the correction 
is purely imaginary. 

As in the NVPT approach based on the complete- 
adiabatic form of the electron-nuclear wave func¬ 
tion M, the first-order perturbation to the electronic 
wave function represents the effect of the non-adiabatic 
coupling between the ground and the excited electronic 
states. Within a fully non-adiabatic approach l[59l[68l . 
it would be possible to compute Eq. [32] . However, it 
has been shown in Ref. ff30l that within DFPT the per¬ 
turbation can be determined by the knowledge of only 
ground-state properties. Eq. (29] is solved by insert¬ 
ing the chosen expression for the electronic wave func¬ 
tion d3Tb and by solving for each order in the perturba¬ 
tion ATfR, t). At the zero-th order we obtain 


Rbo-4°o( r ) ^rM = 0 


( 0 )/ 


(34) 


and at the first-order 


Hbo - 4°o( r ) 


(i) 


V^R, 


^>R( r ) Vv,a. (35) 


Eq. [34] is simply the eigenvalue problem associated to 
the BO Hamiltonian; Eq. [35] is solved in the framework 
of DFPT as illustrated in me Section ll V Bl 

The TDPES of the theory based on the exact factor¬ 
ization remains unaffected if compared to the BO case, 
up to within the first-order perturbation, as shown in 
Eq. [33] . The vector potential, that is identically zero in 
the adiabatic treatment, becomes 


A„(R, t) = —2h Y, Kx (R’ ^) (V 


(o) 


/Pr 




(36) 


This expression is obtained by using Eq. |3l) in the def¬ 
inition of the vector potential given in Eq. [20] . Using 
Eq. [35] in Eq. [36], an alternative expression is derived. 
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that is used in actual calculations, namely 


^(r,*) = -E a ^(R>*)^'( r )' 

v',p 


(37) 


= -2^A£'(R,f) U^p 


v',0 


Hbo ~ eso(R) 


V^R.la 


(38) 


where we have introduced the definition of the matrix 
A%p (R) and the symbol A v n stands for the a Cartesian 
coordinate of the vector potential corresponding to the 
j'-th nucleus. It is instructive to give an alternative for¬ 
mula for the evaluation of the vector potential matrix in 
Eq. [37) , namely 


A& ( R ) = 

e^O 


^eO ,v' (3 (R)c?eO ,i sa (P^) 

W e o(R) 


(39) 


that is obtained by using Eq. |32| in Eq. p8) and acting 
with the BO Hamiltonian on its eigenstates. This expres¬ 
sion is useful to determine the vector potential by com¬ 
bining the NVPT with (explicit) non-adiabatic calcula¬ 
tions. In general, evaluating the vector potential from 
the full electronic wave function in Eq. (llj is difficult 
because the exact electronic state is not known, thus ap¬ 
proximations have to be invoked. Here, we have de¬ 
rived an expression that can instead be used in actual 
calculations. However, in the present paper we focus on 
Eq. [38) and we estimate it within DFPT. 


IV. OBSERVABLES 

A. Current and magnetic dipole moment 

In a time-dependent picture, the expectation values of 
the current and of the magnetic dipole moment on the 
instantaneous state of the system are employed to eval¬ 
uate the rotational strength giving access to the VCD 
spectrum in the linear response regime. We will derive 
their expressions employing the factorized form of the 
full wave function when calculating explicitly the expec¬ 
tation values. 

The current and magnetic dipole moment operators 
are defined as 


N e 


fl = fi e +fi n 


i=1 


= -E^ + EirP- ( 4 °) 


N n 






and 


m 


m 


rn = — 


N e 


N n 


' 2 me 

i= 1 


-r i x 




V =1 


Z v e 
2 M„c 


R„ x P„ 


(41) 

respectively. Here, e is the electronic charge, Z v e is the 
nuclear charge, m and M v are the electronic and nuclear 


masses and c is the speed of light. The position and mo¬ 
mentum operators for the electronic subsystem are indi¬ 
cated as r, and p,, respectively, and similar symbols are 
used for the nuclear operators, R„ and P,,. As expected, 
the vector potential does not appear in Eqs. [40) and [4l) 
since we are not yet calculating an expectation value. 
However, since the nuclear momentum operator in po¬ 
sition representation acts as a derivative w.r.t. the nu¬ 
clear coordinates R, the vector potential appears (only) 
when the derivative acts on the parametric dependence 
of the electronic wave function. Indeed, if the factoriza¬ 
tion is not introduced, such vector potential will never 
be present. 

The expectation values of the operators in Eqs. d40b 
and [IT) on T(r. IE. t) are indicated with the symbol 



= / dR X *(R,f) ($ R (t) U e $a(f) 


N n 


v=l 


(42) 


and 


(m)* = / dRx*(FM) 


($ R (f) |m e |$ R (t)) r 


N, 


m 




x(R,£)- 


(43) 


We will now introduce the following symbols for the 
expectation values of the electronic contributions to the 
current and magnetic dipole moment on the (exact) elec¬ 
tronic wave function. 


A R (t) 

m k(t) 


c l > R (t) b' 3> R (i) 


(4> R (f) |m e |$ R (t)) r 


(44) 

(45) 


If the BO electronic wave function is used to approxi¬ 
mate $ R (r, i), both equations, i.e. the electronic con¬ 
tributions to the expectation values, vanish, as well as 
the vector potential in Eqs. [42) and ( [43) , as mentioned 
above. It is, however, now possible to insert the NVPT 
approximation to the electronic wave function, Eq. [3l) , 
and this leads to the following expressions for the ex¬ 
pectation values. 




iv» 

E /j v e 




P i/ + (R, t) 


(46) 


(m)* 



Pi/ + A„(R, t ) 



(47) 
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Here, we have written the expectation values (on the left 
hand sides) on 'I', the full electron-nuclear wave func¬ 
tion, in terms of expectation values of new observables 
on x, the nuclear wave function only. Therefore, the vec¬ 
tor potential naturally appears in the equations. In ad¬ 
dition, since the electronic wave function has been ap¬ 
proximated, as stated above, by using Eq. ( |3T| ), we ob¬ 
tain that the current and the magnetic dipole moment 
contain an electronic contribution that is first-order (1) 
in the perturbation. The second terms in both equations, 
containing the vector potential, corrects the nuclear con¬ 
tribution to both expectation values and these correc¬ 
tions shall be considered within NVPT since they are 
first-order in the perturbation parameter A„ (R, t) (see 
Eq. {37}). Standard approaches do not consider these 
correction terms, because the vector potential is a quan¬ 
tity that has been introduced only in the context of the 
exact factorization. We will compute explicitly these cor¬ 
rections in Section [V but we can already anticipate that 
while the first (standard) term is 0{ A„), because of the 
P v /M v term, the correction is O ( \„ /M v ) since the vec¬ 
tor potential itself has a linear dependence on the per¬ 
turbation parameter. 

It is worth mentioning here that the advantage of in¬ 
troducing expectation values on the nuclear wave func¬ 
tion, rather than on the full wave function, becomes 
clear when the classical approximation for the nuclear 
subsystem is considered. In this case, due to the prop¬ 
erties of the nuclear wave function in the factorization 
framework (\ is a proper wave function, as it evolves ac¬ 
cording to a TDSE, and leads to the density and current- 
density calculated from the full wave function), the clas¬ 
sical limit can be performed by imposing that the nu¬ 
clear density infinitely localizes, at each time, at the 
classical position denoted by the trajectory. The second 
terms on the RHS of Eqs. 1461 and {47} then become sim¬ 
ply functions of phase-space variables. It is important to 
notice, however, that the vector potential has to be taken 
into account to appropriately relate the nuclear velocity 
and momentum. 


B. Rotational strengths from density functional 
perturbation theory 

The direct numerical solution of Eqs. {34} and {35} is 
very expansive for systems with more than a few de¬ 
grees of freedom. Already the calculation for small chi¬ 
ral molecules requires an approximate treatment of the 
electronic structure problem. In our implementation we 
resort to standard Kohn-Sham (KS) DFT 869M7T1 with 
generalized gradient approximation to the exchange- 
correlation functional 172. 1731 . For simplicity, we will 
limit our discussion to the case of spin saturated closed 
shell systems and drop the explicit notation of the para¬ 
metric dependence on the nuclear positions. 

In the framework of single determinant KS-DFT, 
Eq. {34} directly translates to the standard BO ground- 


state electronic structure problem 


fr(°) (0) 

n KS e o 


^ 0) ( r ) = 0 


(48) 


with the unperturbed KS Hamiltonian H^ s and the un¬ 
perturbed KS orbitals and KS energies ei°' 1 of the 
occupied electronic states o. In DFPT 117414781 , the calcu¬ 
lation of the non-adiabatic correction to the ground-state 
orbitals can be done without explicit knowledge of the 
unoccupied states via a Sternheimer equation |79| 


-P« 


H 


(o) _ ho) 


KS 


£4 1} (r 


P e H^ s [{cf> 0 }]^\ r) (49) 


with a projector on the manifold of the unoccupied 
states P e = 1 — Yh 0 \4 > o)(4’o\ ■ The perturbation Hamil¬ 
tonian on the RHS H^g[{<p 0 }} can depend on the elec¬ 
tronic density response and hence implicitly on the per¬ 
turbed orbitals on the left hand side. This is the case for 
electric field or nuclear displacement perturbations and 
requires a self-consistent solution. Explicitly, Eq. 1491 for 
a nuclear displacement perturbation j reads 


~Pe. 


rr(0) _ (0) 

n KS e o 


p d^\v) _pdHK 1 (0) 

e dRj ~ e ORj 


(50) 


The perturbed KS orbitals da :j ci^ ' 1 (r) are the gradient of 

the KS orbitals </>o( r) w.r.t. a nuclear displacement j. 
They can be used for the calculation of the electronic 
APT in the position form 1131 1301 . 

The corresponding translation of Eq. {35} to DFPT 
reads 


Pe 


H 


(o) 

KS 



PehdR^oX r) Vj. 


(51) 


Also this equation is reminiscent of a Sternheimer equa¬ 
tion. However, instead of an explicit perturbation 
Hamiltonian acting on the unperturbed KS orbitals, the 
RHS is proportional to the gradient of the ground-state 
wave function w.r.t. a nuclear displacement. As al¬ 
ready discussed, this gradient is accessible via Eq. {50} . 
This method requires two response calculations, a self- 
consistent one for the nuclear displacement perturba¬ 
tion and another for the nuclear velocity perturbation. 

Recently, a related approach to the calculation of 
NVPT has been reported l29l which relies on an itera¬ 
tive finite-differences scheme for the construction of the 
intermediate nuclear gradient information. 

With the imaginary correction to the BO electronic 
wave function in Eq. {31} , it is possible to calculate the 
electronic APT £ in the velocity form 


-'a/3 


d{h) 


dR” 




(52) 
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A. (S)-d2-oxirane vs. Oxirane 


and the electronic AAT X 


L a0 


djrhp) 

dK 


2^(0 o |m^|^ a) ) 


(53) 


For the calculation of the magnetic moment, a choice 
of the origin of the position operator has to be made. 
This poses additional complications for the calculation 
of observables in the condensed phase where periodic 
boundary conditions are used. For a detailed discus¬ 
sion, we refer to the literature 1311301 and to Appendix|A| 
The nuclear AAT J is decomposed into its "conven¬ 
tional" contribution and the correction due to the pres¬ 
ence of the vector potential 

Jap = ^ Q/37 f?" + A JZp (54) 


where we used Einstein's summation convention for re¬ 
peated indices. The correction due to the additional 
term in the nuclear magnetic moment in Eq. |47} is given 
by the derivative of 


{Am ^ = 2 Wc e ^ R " AVs = ( 55 ) 


2M„ 


w.r.t. R v a . Written in this form, the correction to the mag¬ 
netic moment depends linearly on the nuclear velocities, 
via the identification A" = /?" . Flowever, this depen¬ 
dence can be removed in the picture of the nuclear AAT. 
To see this, we evaluate the vector potential matrix of 
Eq. <[38j as 

< = 2 E<<(U)I A ks 4 0) I^(E)) (56) 


and take the derivative of Eq. |55| w.r.t. a nuclear veloc¬ 
ity. This gives the correction to the nuclear AAT as 


A ^/3 = 5 ^ 


2 tf,/c 


(57) 


This expression illustrates two features of the correction. 
First, it is non-local in the nuclear contributions, i.e. all 
nuclei contribute to the AAT of a single nucleus. Second, 
the pre-factor contains the inverse nuclear mass, while 
the conventional contribution does not. 


V. NUMERICAL RESULTS 

The presented NVPT has been implemented in our 
development version of the CPMD | [30l [80| electronic 
structure package. The calculations have been per¬ 
formed using DFPT I176U781 with Troullier-Martins lISTI 
pseudo-potentials and the BLYP [72] [73] functional. We 
have employed a plane wave cutoff of 100 Ry. The flu¬ 
orine pseudo-potential with a radius r c = 1.2 has been 
used. The geometry optimizations, harmonic analysis 
and magnetic field perturbation 1201 calculations were 
done using the electronic structure program Gaussian 09 
Revision D.01 H82l employing aug-cc-pVTZ basis set 183] 
and BLYP functional. 


The vector potential from Eq. < [37] has been calculated 
for a small rigid chiral mo 1 ecu 1 e, (S)-d 2 -oxi ra ne shown 
in Fig. [T] As will be clear from the numerical results. 


O 



FIG. 1. (S)-d 2 -oxirane. 


the vector potential contributes only a small fraction to 
the rotational strengths Rk (with k = 1,..., 15 for the 
(S)-d 2 -oxirane and oxirane), as it is computed within a 
perturbation theory approach. The vector potential is 
first-order in the perturbation parameter A„ (R, t) and it 
appears as an explicit ©(M” 1 ) term in the expressions of 
the current and of the magnetic dipole moment. Further 
analysis, currently under investigation, is focussing on 
the calculation of corrections due to the vector potential 
in explicit non-adiabatic molecular dynamics, in order 
to estimate the actual effect of the vector potential on 
observable properties as the VCD signal. 

Before presenting the results for (S)-d 2 -oxirane, let us 
first discuss the case of oxirane, a non-chiral molecule. 
Oxirane differs from (S)-d 2 -oxirane in the deuterium 
atoms, that are replaced by hydrogen atoms. In Fig. [2] 
we draw as blue arrows [54] the velocities correspond^ 
ing to normal modes at 1127 crrU 1 (upper panel) and at 
1489 cm -1 (lower panel), which have been selected as 
examples among the 15 total modes. Perturbations par¬ 
allel to these velocities are used in Eq. d37b to construct 
the vector potential, which is shown as red arrows in the 
figure. It is very interesting to notice that in the case of 




FIG. 2. Vibrational modes at 1127 cm -1 (upper panel) and at 
1489 cm -1 (lower panel) for oxirane, with nuclear velocities 
indicated as blue arrows. The corresponding vector potential 
is shown as red arrows. 

a non-chiral system the vector potential maintains the 
same symmetry of the vibrational modes and is nearly 
anti-parallel to the nuclear displacement: this is what 


8 









one would expect, if the vector potential is not to affect 
the VCD properties, i.e. current and magnetic dipole, 
and thus the rotational strength, of the molecule. 

In the case of (S)-d 2 -oxirane, the results are quite dif¬ 
ferent, as shown in Fig. [3] Again, the velocities cor¬ 
responding to the normalmodes are indicated as blue 
arrows, whereas the vector potential is drawn in red. 
The selected modes are at 896 cm -1 and at 1089 cm” 1 . 
It is clear in this case that (i) a well-defined symmetry 
of the vector potential cannot be identified and, as a 
consequence, (ii) it is not simply (anti-)parallel to the 
normal modes velocities, as was the case for oxirane. 
This behavior thus results in an actual contribution of 
the vector potential to the VCD properties of (S)-d 2 - 
oxirane. Such contribution is quantitatively estimated 




FIG. 3. Vibrational modes at 896 cm -1 (upper panel) and at 
1089 cm -1 (lower panel) for (S)-d 2 -oxirane, with nuclear ve¬ 
locities indicated as blue arrows. The corresponding vector 
potential is shown as red arrows. 


TABLE I. Normal modes for (S)-d 2 -oxirane. The frequencies 
of the modes are indicated in the first column, the rotational 
strengths R are listed in the second column, from Eq. (Rl, the 
corrections A R due to the vector potential are reportedm the 
third (absolute value) and fourth (relative correction) columns. 


V 

(cm” 1 ) (10” 

R 

14 esu 2 cm 2 ) (10” 

A R A R/R 

44 esu 2 cm 2 ) (%) 

647.50 

-0.45 

-0.003 

0.67 

733.42 

10.54 

0.016 

0.15 

769.76 

3.29 

0.001 

0.05 

856.38 

2.70 

0.002 

0.09 

894.67 

-3.89 

0.006 

0.15 

936.33 

-20.26 

0.001 

0.01 

1088.21 

8.34 

-0.027 

0.32 

1093.95 

-4.97 

0.004 

0.09 

1210.44 

10.45 

-0.029 

0.28 

1326.86 

-0.76 

0.0002 

0.03 

1377.38 

-8.17 

0.025 

0.31 

2235.16 

-22.90 

-0.010 

0.04 

2244.19 

16.78 

0.011 

0.07 

3047.68 

-32.59 

-0.063 

0.19 

3054.15 

47.04 

0.047 

0.10 


lows for the calculation of the vector potential beyond 
the NVPT. 


B. Comparison with magnetic field perturbation theory 

Further molecular systems have been investigated, 
namely (R)-propylene-oxide and (R)-fluoro-oxirane 
shown in Figs. |4]an d@ 


by calculating the correction the the rotational strengths 
in Eq. {?} of (S)-d 2 -oxirane, due to the vector potential 
terms in Eqs. {46} and (j47j). Table [I] lists, for all modes in 
the (S)-d 2 -oxirane, these rotational strengths lit and the 
corrections A Rk due to the presence of the vector poten¬ 
tial in the current and in the magnetic dipole moment. 

As discussed above, we notice from the results re¬ 
ported in Table [I] that, despite the fact that the vector 
potential is non-zero, its effect is quite small, being of 
the order In fact, while the M~ x dependence 

in Eqs. {46} and 1471 is removed in the first contribu¬ 
tions, being these first terms proportional to the momen¬ 
tum, the second terms are actually 0(M~ 1 ). We recall, 
however, that in the procedure developed in this pa¬ 
per, the vector potential is evaluated within the NVPT, 
thus being first-order in the perturbation. In a situation 
where the electronic wave function has a strong non- 
adiabatic character, namely where the correction to a 
BO-type wave function is not small in the nuclear veloc¬ 
ity, a larger contribution may be expected. Moreover, in 
the cases where the vector potential is singular, e.g. for 
adiabatic states that are locally degenerate in R-space, 
this correction may become very important. Flowever, 
further studies are required to develop a scheme that al- 


O 



FIG. 4. (R)-propylene-oxide. 



F 


FIG. 5. (R)-fluoro-oxirane. 


In this section we report the dipole D and rotational 
R strengths calculated by employing NVPT, indicated 
with the symbols Dnvp and R^v\> in Tables [II III and |IV| 
and we compare these results with the magnetic field 
perturbation (MFP) theory EOll . Dmfp and Rmfp in the 
tables. Such comparison has been carried out also for 
(S)-d 2 -oxirane (Table [II}. Furthermore, Tables |III| and |TV| 
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TABLE II. Normal modes, dipole and rotational strengths, for 
(S)-d 2 -oxirane. 


z> 

(cm” 1 ) 

DmFP DtsSY p 

(10 4O esu 2 cm 2 ) 

Rmfp -Rnvp 
(10 44 esu 2 cm 2 ) 

647.50 

0.55 

0.85 

-0.35 

-0.45 

733.42 

123.35 

124.88 

8.73 

10.54 

769.76 

53.44 

51.77 

3.17 

3.29 

856.38 

145.31 

145.55 

4.31 

2.70 

894.67 

9.78 

10.24 

-3.37 

-3.89 

936.33 

39.73 

39.24 

-19.14 

-20.26 

1088.21 

3.79 

4.44 

6.95 

8.34 

1093.95 

1.41 

1.71 

-3.98 

-4.97 

1210.44 

26.26 

26.09 

9.56 

10.45 

1326.86 

0.34 

0.37 

-0.91 

-0.76 

1377.38 

11.65 

10.78 

-7.50 

-8.17 

2235.16 

49.17 

50.88 

-22.60 

-22.90 

2244.19 

12.63 

12.81 

16.80 

16.78 

3047.68 

11.43 

11.66 

-32.80 

-32.59 

3054.15 

58.64 

60.16 

46.63 

47.04 


show the corrections AR to the rotational strengths due 
to the vector potential term in Eq. < [3 7) , as alrea dy p re- 
sented for the case of (S)-d 2 -oxirane in Section |V A| In 
all tables the first column indicates the normal modes 
frequency, the second and third columns are the dipole 
strengths from MFP and NVP theories, the forth and 
fifth columns show the rotational strengths from MFP 
and NVP theories. In Tables [ill] and [|V| the sixth and 
seventh columns are the corrections computed from 
Eq. < [37) , which in general are the same order of mag¬ 
nitudeas the corrections reported in Table [I] for (S)-d 2 - 
oxirane. 

From the comparison between the two perturbation 
approaches, we notice an overall very good agreement 
not only in the absolute values of the dipole and rota¬ 
tional strengths, but also in the signs of the rotational 
strengths for the three systems investigated here. The 
MFP theory of Stephens l20l can be considered a "more 
standard" approach, nowadays implemented in most 
quantum-chemistry packages, thus it represents a suit¬ 
able benchmark for the new approach introduced in 
Ref. |30l and discussed in the present work. 


VI. CONCLUSIONS 

One of the main goal of the paper has been to pro¬ 
vide rigorous basis for the development of the NVPT 
approach to VCD. In this context, the complete adiabatic 
approach proposed by Nafie lH9l was adopted in pre¬ 
vious study [30| as starting point, where the electron- 
nuclear wave function is approximated as a single prod¬ 
uct of a (nuclear) vibrational contribution and an elec¬ 
tronic term. In particular, such electronic term contains 
corrections to the BO state which are first-order in the 
nuclear velocity. In the present work, we make this idea 


TABLE III. Normal modes, dipole and rotational strengths 
(with corrections), for (R)-propylene-oxide. AR and AR/R in 
the last columns are indicated in 10 _44 esu 2 cm 2 and %, respec¬ 
tively. 


(cm” 1 ) 

-Dmfp Dnvp 
(10 40 esu 2 cm 2 ) 

.Rmfp Rnvp 
(10 44 esu 2 cm 2 ) 

AR A R/R 

202.12 

6.91 

7.12 

3.54 

3.47 

-0.001 

0.03 

355.33 

45.22 

46.63 

-12.84 

-12.56 

0.008 

0.06 

398.40 

38.74 

38.86 

-3.72 

-3.79 

0.004 

0.11 

717.36 

55.11 

52.51 

13.88 

13.21 

0.005 

0.04 

795.26 

217.23 

219.15 

2.47 

1.62 

0.007 

0.44 

875.72 

18.54 

17.25 

26.36 

26.99 

0.039 

0.14 

929.21 

51.55 

51.53 

-35.52 

-37.03 

-0.049 

0.13 

1008.29 

24.84 

26.62 

2.88 

4.53 

-0.004 

0.09 

1089.09 

18.53 

19.17 

-6.03 

-6.56 

0.006 

0.09 

1112.88 

7.92 

7.79 

6.65 

7.50 

0.023 

0.31 

1126.68 

11.68 

12.56 

-13.44 

-14.67 

-0.034 

0.23 

1150.27 

1.51 

1.40 

1.54 

1.23 

0.003 

0.24 

1246.96 

19.77 

19.85 

-8.06 

-8.01 

-0.004 

0.05 

1371.08 

10.35 

9.83 

3.30 

3.53 

0.007 

0.19 

1388.57 

60.08 

60.10 

13.99 

15.15 

0.007 

0.05 

1447.69 

13.15 

14.16 

1.34 

1.45 

0.005 

0.32 

1461.62 

15.41 

16.62 

-1.69 

-1.90 

-0.008 

0.42 

1480.79 

10.14 

9.99 

4.66 

4.69 

-0.005 

0.11 

2955.51 

27.68 

28.86 

1.64 

1.64 

0.0002 

0.01 

3000.54 

41.29 

44.50 

-0.29 

0.20 

-0.009 

4.53 

3005.59 

22.70 

24.14 

5.13 

6.04 

-0.034 

0.56 

3007.28 

22.57 

22.86 

-13.92 

-15.14 

0.053 

0.35 

3032.47 

47.06 

50.07 

7.29 

7.16 

-0.019 

0.27 

3079.38 

41.31 

41.09 

-7.19 

-7.31 

0.013 

0.17 


exact, in the sense that the starting point is not an ap¬ 
proximate factorized form of the full wave function. The 
starting point is provided by the exact factorization of 
the electron-nuclear wave function, where approxima¬ 
tions are inserted at a later stage in order to make nu¬ 
merical calculations feasible. The method outlined here 
can thus be seen as a rigorous basis for NVPT: at the first 
stage of the derivation we describe how to recover the 
BO working equation from the exact electronic equation 
and at the second stage a perturbation to BO is consid¬ 
ered. Also, this perturbation does not rely on the use 
of the nuclear velocity as small parameter, in fact such 
parameter is, more generally, related to the spatial vari¬ 
ations of the nuclear wave function from the factoriza¬ 
tion. Only in the classical limit, at O(h 0 ), these varia¬ 
tions lead to an interpretation in terms of nuclear veloc¬ 
ity. In the new approach presented here, a full quantum 
picture can be maintained, without invoking the classi¬ 
cal approximation. 

The second main result confirms the importance of 
using the exact factorization as starting point for the 
development of approximations. The time-dependent 
vector potential of the theory naturally appears in the 
observables, i.e. the current and the magnetic dipole 
moment, necessary for the calculation of the VCD spec¬ 
trum. Therefore, within the perturbation approach pre- 
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TABLE IV. Normal modes, dipole and rotational strengths 
(with corrections), for (R)-fluoro-oxirane. A R and A R/R in 
the last columns are indicated in 10 _ 44 esu 2 cm 2 and %, respec¬ 
tively. 


z> 

(cm- 1 ) 

Dmfp T>nvp 
(10 40 esu 2 cm 2 ) 

.Rmfp Rnvp 
(10 44 esu 2 cm 2 ) 

A R AR/R 

411.61 

52.63 

53.11 

9.48 

9.80 

-0.003 

0.03 

482.91 

30.46 

31.23 

-3.10 

-2.91 

0.002 

0.06 

733.56 

124.68 

123.64 

40.79 

39.91 

0.031 

0.08 

804.61 

501.12 

497.82 

-12.79 

-9.85 

0.007 

0.07 

927.57 

244.98 

246.52 

-27.57 

-34.46 

-0.052 

0.15 

1059.05 

28.75 

25.77 

-9.38 

-9.09 

-0.019 

0.21 

1069.68 

312.66 

321.09 

22.55 

22.47 

0.048 

0.21 

1106.47 

5.52 

4.95 

-8.37 

-8.84 

-0.022 

0.25 

1125.52 

11.52 

11.28 

4.11 

4.77 

0.005 

0.10 

1252.78 

88.68 

87.54 

-0.07 

1.73 

0.006 

0.32 

1344.65 

150.66 

150.18 

-6.39 

-7.04 

-0.016 

0.23 

1470.12 

42.55 

44.05 

0.73 

0.77 

0.008 

1.06 

3024.87 

20.22 

21.53 

1.64 

1.53 

0.003 

0.16 

3068.24 

22.58 

23.21 

-1.07 

-1.00 

0.008 

0.84 

3115.60 

14.50 

14.02 

0.34 

0.37 

-0.007 

1.79 


sented in the paper, we have evaluated the vector po¬ 
tential using the harmonic approximation for the nu¬ 
clear motion. In this case, the contribution has been 
shown to be small, but only further investigation, for 
instance in the context of non-adiabatic molecular dy¬ 
namics, will clarify the actual extent of non-adiabatic 
corrections to the VCD signal. Also, situations where 
the non-adiabatic couplings are important shall be in¬ 
vestigated, for instance for low-lying excited states l85l , 
where the exact factorization approach offers a strategy 
to overcome the limitations of BO approximation in a 
rigorous way. 

According to the procedure presented in this work, 
NVPT is suitable for an implementation in any ab-initio 
molecular dynamics code. Therefore, NVPT can be eas¬ 
ily employed for the study of VCD properties of chiral 
molecules in solutions and for direct comparison with 
experimental data. Such procedure allows also to eval¬ 
uate the corrections due the vector potential from the 
exact factorization approach. As it requires a DFPT cal¬ 
culation for each geometry sampled by the molecular 
dynamics trajectory, the numerical cost of a NVPT cal¬ 
culation is slightly larger than standard BO molecular 
dynamics but indeed feasible. 
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Appendix A: Invariance under choice of the origin 

One of the main problems connected to the evaluation 
of molecular properties and spectroscopies depending 
on the magnetic field is to assure origin invariance of 
the final results. In case of VCD, this requires the eval¬ 
uation of the electric and magnetic dipole moments, or 
accordingly the APT and the AAT. While the APT shows 
no origin dependency, the exact AAT transforms under 
shifts of the origin 0 = 0' + A as 

M v ° p = Ks- (Al) 

7<5 

The rotational strength as a physical observable is gauge 
invariant 

Rk= E sz k s£ k 

a a'(3 vis' 

E (A2) 

olol'{3 7<5 vv' 

since the second terms constitute triple products con¬ 
taining two identical vectors. 

The evaluation of origin dependent operators under 
periodic boundary conditions has been extensively dis¬ 
cussed in the literature 186H881 . A convenient approach 
is the combination of statewise origins with maximally 
localized Wannier orbtials, which has been applied suc¬ 
cessfully to the calculation of nuclear magnetic reso¬ 
nance chemical shifts l89l [901- The canonical <j> 0 and 
localized ip a states are mutually related via the uni¬ 
tary transformation for the unperturbed ground-state 
orbitals 

Wo) =J2 U oo'\M. (A3) 

o' 

This approach is based on the natural assumption that 
the response orbitals are sufficiently localized in the re¬ 
gion of their respective unperturbed ground-state or¬ 
bitals. In the distributed origin (DO) gauge, the position 
operators are calculated with the corresponding Wan¬ 
nier center as its statewise origin 

v 0 = (<p 0 |f|Po)- (A4) 

The electronic AAT in a statewise DO gauge then is 
given by 

Waff )do = - r oi)PS^-ysWo}^ a) ) (A5) 

and can be translated back to the common origin form 
via 

V$ = 

O 

+ (A6) 

07 6 
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where £'f is the contribution of the state o to the elec- culation are the same for canonical and Wannier or- 
tronic APT. The numerical results in a supercell cal- bitals 1130] . 
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